$\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\br}[1]{\left[#1\right]}$ $\newcommand{\cbr}[1]{\left\{#1\right\}}$ $\newcommand{\abs}[1]{\left|#1\right|}$
Chaos is aperiodic long-term behavior in a deterministic system that exhibits sensitive dependence on initial conditions.
Any study of chaos begins with the Lorenz equations:
\[ \begin{align*} \dot{x} &= \sigma\pa{y-x}, \\ \dot{y} &= x\pa{r-z} - y, \\ \dot{z} &= xy-bz, \end{align*}\]for parameters $\sigma, r, b > 0$: $\sigma$ is the Prandtl number, $r$ is the Rayleigh number, and $b$ is nameless. As Lorenz did, we will study the particular case $\sigma = 10$, $b = 8/3$, and $r = 28$. Note: From here on forward I will be referring to the Lorenz equations as system (1).
In this code, we seek to plot trajectories in three dimensions and create the strange attractor, known as the Lorenz attractor or Lorenz butterfly (which has a fractal dimension between 2 and 3).
Simply from glancing at the equations in sytem (1) we see the following fixed point: $\pa{x^\ast, y^\ast, z^\ast} = \pa{0,0,0}$, which holds $\forall\,\sigma, r, b$. For $r>1$, another fixed point can be quickly discovered:
\begin{gather*} 0 = \sigma\pa{y-x} \Rightarrow y^\ast = x^\ast; \\ 0 = x\pa{r-z} - y = x\pa{r-z-1} \Rightarrow\boxed{z^\ast = r-1}; \\ 0 = xy-bz = x^2-bz \Rightarrow\boxed{x^\ast = y^\ast = \pm\sqrt{b\pa{r-1}}}. \end{gather*}We will call these fixed points $C^+$ and $C^-$. As $r\to 1^+$, $C^+$ and $C^-$ converge to the origin to form a pitchfork bifurcation.
The linearization at the origin is
\begin{align*} \dot{x} &= \sigma\pa{y-x}, \\ \dot{y} &= rx - y, \\ \dot{z} &= -bz, \end{align*}obtained by omitting nonlinearities in system (1). Since the third equation is now decoupled, we see that the solution is $z(t) = z_0\exp{\pa{-bt}}$, where $z_0$ is the initial condition. The other degrees of freedom are governed by the system
$$\begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} -\sigma & \sigma \\ r & -1 \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix},$$with $\tau = -\pa{\sigma + 1} < 0$ and $\Delta = \sigma\pa{1-r}$. If $r>1$, $\Delta < 0$ $\Rightarrow$ the origin is a saddle point. If $r<1$, the origin is a sink. Note that $\tau^2 - 4\Delta = \pa{\sigma +1}^2 - 4\sigma\pa{1-r} = \pa{\sigma - 1}^2 + 4\sigma r > 0$, so the origin is a stable node for $r<1$.
To show that the origin is globally stable, we construct a Lyapunov function. Consider $V: \mathbb{R}^3\to\mathbb{R}$, defined by
$$V\pa{x, y, z} \equiv \frac{1}{\sigma}x^2+y^2+z^2.$$Then,
\begin{align*} \dot{V} &= 2\pa{\frac{1}{\sigma}x\dot{x} + y\dot{y} + z\dot{z}}, \\ \frac{1}{2}\dot{V} &= x\pa{y-x} + y\br{x\pa{r-z} - y} + z\pa{xy-bz} \\ &= \pa{1+r}xy - x^2 - y^2 - bz^2 \\ &= -\pa{x-\frac{1+r}{2}y}^2 - \br{1-\pa{\frac{1+r}{2}}^2}y^2 - bz^2, \end{align*}we see that $\dot{V}$ is strictly negative $\forall\, x,y,z$, so long as $r<1$. $\dot{V} = 0$ only when $\pa{x,y,z} = \pa{0,0,0}$, therefore the origin is globally stable for $r<1$.
Suppose $r>1$ so that $C^+$ and $C^-$ exist. We can determine their stability from the eigenvalues of the Jacobian matrix for the system.
$$J = \begin{pmatrix} \partial_x\dot{x} & \partial_y\dot{x} & \partial_z\dot{x} \\ \partial_x\dot{y} & \partial_y\dot{y} & \partial_z\dot{y} \\ \partial_x\dot{z} & \partial_y\dot{z} & \partial_z\dot{z} \end{pmatrix} = \begin{pmatrix} -\sigma & \sigma & 0 \\ r-z & -1 & -x \\ y & x & -b \end{pmatrix} = \begin{pmatrix} -\sigma & \sigma & 0 \\ 1 & -1 & \mp\sqrt{b\pa{r-1}} \\ \pm\sqrt{b\pa{r-1}} & \pm\sqrt{b\pa{r-1}} & -b \end{pmatrix},$$after replacing $x, y, z$ with $x^\ast, y^\ast, z^\ast$, respectively. Then,
\begin{align*} 0 &= \text{det}\pa{J-\lambda I_3} \\ &= \text{det}\begin{pmatrix} -\sigma-\lambda & \sigma & 0 \\ 1 & -1-\lambda & \mp\sqrt{b\pa{r-1}} \\ \pm\sqrt{b\pa{r-1}} & \pm\sqrt{b\pa{r-1}} & -b-\lambda \end{pmatrix} \\ &= \lambda^3 + \pa{\sigma + 1 + b}\lambda^2 + b\pa{\sigma + r}\lambda + 2b\sigma\pa{r - 1}, \end{align*}let $\lambda = i\omega$, so that
\begin{equation} \tag{2} 0 = -i\omega^3 - \pa{\sigma + 1 + b}\omega^2 + b\pa{\sigma + r}i\omega + 2b\sigma\pa{r - 1}. \end{equation}Both the imaginary and real parts must be zero, we see from the imaginary component that
\begin{align*} 0 &= -i\omega^3 + b\pa{\sigma + r}i\omega, \\ \therefore \omega &= 0, ~ \omega^2 = b\pa{\sigma + r}; \end{align*}but $\omega \neq 0$ because then Equation (2) would not hold, so $\omega^2 = b\pa{\sigma + r}$ is the only solution. We substitute this result into the real component
\begin{align*} 0 &= - \pa{\sigma + 1 + b}\omega^2 + 2b\sigma\pa{r - 1} \\ &= - \pa{\sigma + 1 + b}\br{b\pa{\sigma + r}} + 2b\sigma\pa{r - 1} \\ &= -\pa{b + b^2 - b\sigma}r - \pa{b\sigma^2 + b^2\sigma + 3b\sigma}, \\ \therefore r &= \sigma\pa{\frac{\sigma + b + 3}{\sigma - 1 - b}}; \end{align*}note that $r>0 ~ \Rightarrow ~ \sigma > 1 + b$. Let $$r_H \equiv \frac{\sigma\pa{\sigma + b + 3}}{\sigma - 1 - b},$$ where we use the subscript $H$ to denote $C^+$ and $C^-$ losing their stability in a Hopf bifurcation at this value.
import numpy as np # For math operations and variable storing
from scipy.integrate import odeint # To solve the differential equations
We start by defining all of our variables. In particular, we define $\sigma = 10$, $r = 28$, and $b = 8/3$, as mentioned above. We also define the number of steps in our iteration with the steps variable. Our initial conditions are stored in the X0 and X1 variables, and t is our time array.
sigma = 10
r = 28
b = 8/3
steps = 10000
t = np.linspace(0, 100, steps)
X0 = [0.0, 1.0, 1.05]
X1 = [0.0, 1.0, 1.0]
# Compute the fixed points
x_fixed = [0, np.sqrt(b * (r - 1)), -np.sqrt(b * (r - 1))]
y_fixed = [0, np.sqrt(b * (r - 1)), -np.sqrt(b * (r - 1))]
z_fixed = [0, r - 1, r - 1]
We now define a function lorenz which describes our system of equations in group (1). It takes in the initial conditions and all parameters in the Lorenz equation.
def lorenz(X, t, sigma, r, b):
x, y, z = X
dxdt = sigma * (y - x) # First Lorenz equation
dydt = x * (r - z) - y # Second Lorenz equation
dzdt = x * y - b * z # Third Lorenz equation
return [dxdt, dydt, dzdt]
To solve the differential equations we use odeint. This decision was made in order to improve accuracy and reduce computational complexity. From the figures that I've generated, I prefer the ones produced by odeint over standard numerical techniques like Forward Euler.
solution0 = odeint(lorenz, X0, t, args=(sigma, r, b))
# Extracting components from the solution
x0 = solution0[:, 0]
y0 = solution0[:, 1]
z0 = solution0[:, 2]
solution1 = odeint(lorenz, X1, t, args=(sigma, r, b))
# Extracting components from the solution
x1 = solution1[:, 0]
y1 = solution1[:, 1]
z1 = solution1[:, 2]
In this section we plot some amazing figures. Note that most of the code is commented out to decrease the size of this file.
The following cell imports all of the necessary packages to plot using plotly.
import style # Custom styling for figures
import plotly.graph_objects as go # To create the animations
import plotly.io as pio
pio.templates["custom_template"] = style.template
pio.templates.default = "custom_template"
This cell sets the automatic perspective to face the front of the Lorenz attractor.
camera = dict(
eye=dict(x=1, y=-1, z=0.2)
)
The following cell prints the Lorenz attractor for a single data set.
# Define labels for fixed points
labels_fixed = ['Origin', 'Fixed Point 1', 'Fixed Point 2']
fig = go.Figure(data=[
# Plotting the trajectory
go.Scatter3d(
x=x0, y=y0, z=z0, mode='lines', line=dict(width=0.5), name='Trajectory'
),
# Plotting the fixed points
go.Scatter3d(
x=x_fixed, y=y_fixed, z=z_fixed, mode='markers', marker=dict(size=5, color="#9FFFCB"),
text=labels_fixed, hoverinfo='text', name='Fixed Points'
)
])
fig.update_layout(
title=f"Lorenz Attractor for sigma = {sigma}, r = {r}, and b = {b:.2g}",
scene=dict(
xaxis_title='x axis',
yaxis_title='y axis',
zaxis_title='z axis',
camera=camera
))
fig.show() # To print the figure in the notebook
Strogatz, Steven. "Lorenz Equations." Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. CRC Press, 2018. pp. 309 - 328.